#ifndef BIVARIATEP_HPP
#define BIVARIATEP_HPP

#define MAXITER 10

class Bivariatep{
	public:
		Bivariatep(double thres)
		{
			_thres=thres;
			d=new boost::math::normal_distribution<>(0,1); 

		}
		double operator()(double rho, double b1, double b2)
		{
			double res=(cdf(*d,b1)+cdf(*d,b2))/2;
			res+=E(b1,b2-b1*rho/(b1*sqrt(1-rho*rho)));
			res+=E(b2,b1-b2*rho/(b2*sqrt(1-rho*rho)));
			if(b1*b2<0 |((b1*b2!=0) | (b1+b2 < 0))) res+=0.5;
				
			return res;
		}
		double E(double b, double bh)
		{
			double prec=10;
			double curr=atan(bh);
			double a=0;
			double c=0;
			double j=0;
			int fact=1;
			double sum=0;
			while(abs(prec-curr)>_thres | j> MAXITER)
			{
				prec=curr;
				if(j!=0) fact*=j;
				sum+=pow(b,2*j)/pow(2,j)*fact;	
				if((double)j/2==(int)(j/2))
				{
					c=(double)1/(2*j+1)*(1-exp(b*b/2)*sum);
				}else{
					c=-(double)1/(2*j+1)*(1-exp(b*b/2)*sum);
				}
				curr-=c*pow(bh,2*j+1);
				cout <<"curr " << curr;
				j++;
			}
			return curr/(2*PI);

		}
	private:
		boost::math::normal_distribution<> *d;
		double 	_thres;

};
#endif
